%无量纲化万有引力常量
%G = 1
%地月质量和6.04586*10^24Kg = 1
%地月轨道半径384400Km = 1
%地月环绕平均速度1.024544Km/s = 1
%特征时间=特征长度/特征速度=375191.305s=4.34249195day = 1

clear
clc
close all
h=1e-5;%步长，可以调整精确程度，先大后小
t=0:h:5;%时间
global m
m = 76 * pi / 180;%飞船开始时地月与长轴的角度
n = -45 * pi / 180;%飞船离开地球时的角度
e = 0.004567;%月球自转速度
u = 0.452213;%地球自转速度大小
v = 10.38; %飞船自带初速度大小
theta = 45 * pi / 180; %飞船自带初速度方向
mxyuv1=[0.987849 -0.012151*cos(m) -0.012151*sin(m) 0.012151*sin(m) -0.012151*cos(m)]';%地球
mxyuv2=[0.012151 0.987849*cos(m) 0.987849*sin(m) -0.987849*sin(m) 0.987849*cos(m)]';%月球cos(m) sin(m)
mxyuv3=[0.165402441*10^-24 0.016573881*cos(n)-0.012151*cos(m) 0.016573881*sin(n)-0.012151*sin(m) -u*sin(n)+v*cos(theta)+0.012151*sin(m) u*cos(n)+v*sin(theta)-0.012151*cos(m)]';%飞船，质量只需要满足数量级


%计算双星系统，定性计算，没有代入实际物理参数
figure('name','飞船位移','position',[100,200,500,500])
global h1
h1 = animatedline('color','b');
global h2
h2 = animatedline('color','r');
global h3
h3 = animatedline('color','y');%黄色是飞船，蓝色是地球，红色是月球
colordef black
grid on
axis equal
view(2)
xlabel('X');
ylabel('Y');
hold on

figure('name','能量','position',[600,200,500,500])
global h4
h4 = animatedline('color','g');%绿色，动能
global h5
h5 = animatedline('color','m');%品红，势能
global h6
h6 = animatedline('color','y');%黄色，总能量
grid on
axis equal
view(2)
xlabel('X');
ylabel('Y');
hold on

figure('name','速度分量','position',[1100,200,500,500])
global h7
h7 = animatedline('color','c');%x方向，亮青
global h8
h8 = animatedline('color','m');%y方向，品红
global h9
h9 = animatedline('color','y');%黄色，总能量
grid on
axis equal
view(2)
xlabel('X');
ylabel('Y');
hold on


y=ODE_RK4_hyh(t,h,[mxyuv1;mxyuv2;mxyuv3]);

%plot(y(2,:),y(3,:),y(7,:),y(8,:))%绘轨迹图

function y=ODE_RK4_hyh(x,h,y0)
%4阶RK方法
%h间隔为常数的算法
global m
global h1
global h2
global h3
global h4
global h5
global h6
global h7
global h8
global h9
y=zeros(size(y0,1),size(x,2));
y(:,1)=y0;
for ii=1:length(x)-1
    yn=y(:,ii);
    xn=x(ii);
    K1=Fdydx(xn,yn);
    K2=Fdydx(xn+h/2,yn+h/2*K1);
    K3=Fdydx(xn+h/2,yn+h/2*K2);
    K4=Fdydx(xn+h,yn+h*K3);
    y(:,ii+1)=yn+h/6*(K1+2*K2+2*K3+K4);%RK4方法
    if ii/1000 == 94.6%第二幅图能量的黄线代表的时间
        y(14,ii+1)=0.362;
        y(15,ii+1)=-1.493;%变轨
    end
    if mod(ii,100) == 50 || ii == 1 || ii == length(x)-1
      addpoints(h1, y(2,ii),y(3,ii));
      addpoints(h2, y(7,ii),y(8,ii));
      addpoints(h3, y(12,ii),y(13,ii));
      addpoints(h4, ii/1000,0.5*(y(14,ii)^2+y(15,ii)^2))
      addpoints(h5, ii/1000,0.987849/(((y(2,ii)-y(12,ii))^2+(y(3,ii)-y(13,ii))^2))^0.5+0.012151/(((y(7,ii)-y(12,ii))^2+(y(8,ii)-y(13,ii))^2))^0.5)
      addpoints(h6, ii/1000,0.5*(y(14,ii)^2+y(15,ii)^2)-0.987849/(((y(2,ii)-y(12,ii))^2+(y(3,ii)-y(13,ii))^2))^0.5-0.012151/(((y(7,ii)-y(12,ii))^2+(y(8,ii)-y(13,ii))^2))^0.5)%加点
      addpoints(h7, ii/1000 ,y(14,ii));
      addpoints(h8, ii/1000 ,y(15,ii));%加点    
      addpoints(h9, ii/1000,0.5*(y(14,ii)^2+y(15,ii)^2)-0.987849/(((y(2,ii)-y(12,ii))^2+(y(3,ii)-y(13,ii))^2))^0.5-0.012151/(((y(7,ii)-y(12,ii))^2+(y(8,ii)-y(13,ii))^2))^0.5)%加点
      drawnow limitrate
        if (y(12,ii)-y(7,ii))^2+(y(13,ii)-y(8,ii))^2 < 0.00451977^2 || (y(12,ii)-y(2,ii))^2+(y(13,ii)-y(3,ii))^2 < 0.01657388^2
            addpoints(h1, y(2,ii),y(3,ii));
      addpoints(h2, y(7,ii),y(8,ii));
      addpoints(h3, y(12,ii),y(13,ii));
      addpoints(h4, ii/1000,0.5*(y(14,ii)^2+y(15,ii)^2))
      addpoints(h5, ii/1000,0.987849/(((y(2,ii)-y(12,ii))^2+(y(3,ii)-y(13,ii))^2))^0.5+0.012151/(((y(7,ii)-y(12,ii))^2+(y(8,ii)-y(13,ii))^2))^0.5)
      addpoints(h6, ii/1000,0.5*(y(14,ii)^2+y(15,ii)^2)-0.987849/(((y(2,ii)-y(12,ii))^2+(y(3,ii)-y(13,ii))^2))^0.5-0.012151/(((y(7,ii)-y(12,ii))^2+(y(8,ii)-y(13,ii))^2))^0.5)
      addpoints(h7, ii/1000 ,y(14,ii));
      addpoints(h8, ii/1000 ,y(15,ii));%加点
      addpoints(h9, ii/1000,0.5*(y(14,ii)^2+y(15,ii)^2)-0.987849/(((y(2,ii)-y(12,ii))^2+(y(3,ii)-y(13,ii))^2))^0.5-0.012151/(((y(7,ii)-y(12,ii))^2+(y(8,ii)-y(13,ii))^2))^0.5)
      drawnow limitrate
      break
        end
    end
end
end

function dydx=Fdydx(x,y)
%将原方程整理为dy=F(y,x)的形式
N=numel(y)/5;%N个球体
G=1;%无量纲化
%计算星球之间的引力
m=y(1:5:end);
x0=y(2:5:end);
y0=y(3:5:end);
[x1,x2]=meshgrid(x0,x0);
[y1,y2]=meshgrid(y0,y0);
[m1,m2]=meshgrid(m,m);%meshgrid，之前讲过，创建二维网格
dx=x1-x2;
dy=y1-y2;
ax=-G.*dx./(dx.^2+dy.^2).^(1.5).*m2;
ay=-G.*dy./(dx.^2+dy.^2).^(1.5).*m2;
for k=1:N
    ax(k,k)=0;
    ay(k,k)=0;%自己对自己不形成引力
end

%建立dydx
dydx=zeros(size(y));
%1 质量不变
dydx(1:5:N*5-4)=0;
%2 x导数
dydx(2:5:N*5-3)=y(4:5:N*5-1);
%3 y导数
dydx(3:5:N*5-2)=y(5:5:N*5-0);
%4 x加速度
dydx(4:5:N*5-1)=sum(ax,1);
%3 y加速度
dydx(5:5:N*5-0)=sum(ay,1);
end
